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I  Introduction 

This  is  the  first  annual  report  of  the  atmospheric  wave 
studies  project  at  the  Stanford  Geophysics  Department.  The  prin- 
clpal  objective  is  to  understand  atmospheric  pressure  fluctuations 
in  the  period  range  1-15  minutes.  At  present  we  cannot  explain 
most  of  the  observed  variation  at  Stanford.  Next  year  when  we 
have  three  Instruments  Installed  it  is  likely  that  we  will  be 
able  to  ascribe  much  of  the  fluctuation  to  jet  stream  instabili¬ 
ties.  This  has  been  the  case  for  data  recorded  in  Boston  (Madden 
and  Claerbout)  and  New  York  (Tolstoy) ,  however  the  jet  stream  is 
much  farther  south  in  the  western  U.  S.  and  we  will  not  pre-judge 
thin  matter. 

During  the  past  year  there  have  been  atmospheric  nuclear 
explosions  set  off  in  China  and  the  South  Pacific  (by  the  French). 

The  South  Pacific  explosions  were  too  weak  and  far  away  to  be 
recorded  on  our  microphone.  These  explosions  offer  a  unique  sit¬ 
uation  for  atmospheric  studies.  The  sound  waves  emitted  from 
such  an  explosion  sample  winds  and  temperatures  ut  altitudes, 
latitudes  and  longitudes  which  are  Inaccessible  to  conventional 
meteorological  observation.  The  theoretical  section  of  this 
report  represents  an  attempt  to  develop  a  means  of  simulating 
the  wave '8  propagation  about  the  earth. 
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II  Antipodal  Focus  on  an  Irregular  Earth 

The  sound  waves  emitted  from  nuclear  explosions,  exploding 
volcanos,  and  great  earthquakes  have  been  observed  to  propagate 
more  than  once  around  the  earth.  Fcr  example,  Wexler  and  Hess 
(1962,  Global  Atmospheric  Pressure  Effects  of  the  October  30,  1961, 
Explosion,  Journal  of  Geophysical  Research,  v.  67,  n.  10,  p.  3875- 
3887)  summarize  data  from  around  the  world  for  an  explosion  In  Si¬ 
beria.  The  disturbance  seen  at  the  antipodal  point  is  of  special 
interest  because  from  an  over  Idealized  point  of  view  one  expects 
that  this  disturbance  could  reach  a  magnitude  comparable  to  the 
magnitude  at  the  source.  Besides  the  obvious  need  to  account  for 
dissipative  effects,  observations  indicate  the  extreme  Importance 
of  horizontal  atmospheric  variations  in  defocuslng  the  antipodal 
focus .  See  the  highly  Irregular  arrival  time  maps  in  figures  1 
to  3.  By  looking  at  the  U.  S.  one  realizes  that  complexity  Increases 
with  the  density  of  observation. 

There  has  been  much  theoretical  work  on  global  air  waves  In  the 
1960's  but  It  seems  that  perhaps  none  of  It  has  addressed  itself 
to  the  horizontal  atmospheric  variations  which  are  so  apparent  In 
these  observations.  Our  objective  in  this  study  Is  to  develop  a 
feasible  mathematical-computational  scheme  for  quantitative  descrip¬ 
tion  of  these  phenomena.  We  would  hope  to  be  able  to  answer  such 
questions  as:  Which  is  more  important,  temperature  variations  or 
wind  variations?  From  observation  of  the  waves  may  we  deduce  any 


Figure  1:  Arrival  time  of  primary  wave  over  the  northern  hemisphere. 


(from  Wexler  and  Hass) 
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Figure  2:  Arrival  time  of  the  primary  wave  over  the  southern  hemisphere 
(from  Wexler  and  Hass) 


Figure  3:  Arrival  time  of  the  wave  returning  to  the  northern  hemisphere 
after  going  through  the  antipodes,  (from  Wexler  and  Hass) 


atmospheric  parameters  which  are  not  observed  by  the  world  wide 
meteorological  network?  If  we  knew  the  meteorological  parameter* 
how  accurately  could  we  estimate  the  source  strength  and  location? 

In  theoretical  work  there  seems  to  be  a  principle  of  compounding 
complexity.  If  one  attempts  to  enlarge  a  theory  to  include  a  more 
complicated  model  then  the  new  complication  interacts  with  all  the 
old  complications  to  create  a  theory  which  is  too  unwieldy  to  attempt 
to  program.  Our  approach  in  this  matter  will  be  to  throw  out  some 
of  the  complexities  of  the  theory  developed  by  Pfeffer  and  Zarlchney, 
Press  and  Harkrlder,  Pierce,  Tolstoy,  Claerbout,  and  others  in  the 
hope  of  being  able  to  include  the  horizontal  atmospheric  variations 
responsible  for  figures  1-3.  With  luck  what  we  throw  out  should  be 
less  important  than  what  we  include  and  the  final  result  should  be 
programmable  with  reasonable  effort.  We  begin  by  throwing  out  the 
altitude  variation  of  temperature  and  wind.  These  will  be  replaced 
by  some  kind  of  average.  The  average  would  probably  be  weighted  ac¬ 
cording  to  the  energy  density  of  the  so-called  Lamb  wave. 

In  a  homogeneous  constant  temperature  atmosphere  this  is  an  expc'.en- 
tlal  with  a  scale  height  of  15  to  30  kilometers.  From  these  assump¬ 
tions  we  may  derive  a  scalar  wave  equation  in  two  dimensions,  that 
is,  on  the  earth's  surface.  At  this  point  we  could  in  concept  solve 
the  equation  by  finite  differences  on  a  grid.  Except  for  the  spher¬ 
ical  geometry  the  actual  calculation  would  not  differ  a  great  deal 
from  the  selsmological  calculations  of  David  M.  Boore  in  his  recent 
Ph.D  thesis  at  MIT.  In  practice  a  15  minute  period  wave  sampled 
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at  10  points  per  wavelength  would  require  a  megaword  computer  memory 
to  represent  the  earth.  Consideration  of  concepts  of  numerical  holo¬ 
graphy  Indicates  that  about  three  orders  of  magnitude  of  computational 
simplification  Is  possible.  The  full  blow.:  two  dimensional  wave 
equation  has  the  capability  of  describing  a  superposition  of  many 
waves  each  traveling  in  an  arbitrary  compass  direction.  As  the 
data  shows  we  are  dealing  with  merely  a  single  outgoing  quasl-clrcular 
wave  front.  The  full  generality  of  the  two  dimensional  wave  equation 
is  not  required  and  we  are  at  liberty  to  seek  a  simpler  controllng 
equation.  From  a  ray  point  of  view  we  may  say  that  the  atmospheric 
Irregularity  is  sufficient  to  cause  rays  to  bend  and  perhaps  cross 
over  one  another  but  any  reflections  set  up  by  the  irregularity  may 
be  neglected.  Thus  the  situation  is  even  simpler  than  reflection 
seismology,  a  field  where  similar  concepts  have  been  applied  (Claer- 
bout,  1970,  Coarse  grid  calculations  of  waves  in  Inhomogeneous  media 
with  application  to  delineation  of  complicated  seismic  structure, 
Geophysics,  v.  35,  n.  2  or  3  and  Claerbout,  1970,  Toward  a  unified 
theory  of  reflector  mapping,  Geophysics,  in  press).  Physically  a 
very  natural  description  is  to  strrt  at  the  source  and  extrapolate 
the  disturbance  out  around  the  earth  as  many  times  around  the  earth 
as  it  can  be  observed.  Mathematically  we  must  consider  a  single 
frequency  at  a  time.  Mathematically  speaking  the  total  pressure 
disturbance  is  periodic  with  distance  from  the  source,  the  period 
being  equal  to  the  earth's  circumference.  What  we  really  need  is 
an  equation  governing  the  out-going  wave,  not  the  scalar  Helmholtz 


equation  which  governs  the  total  pressure  at  a  single  frequency  where 

"total"  means  the  sum  of  all  the  waves  regardless  of  how  many  times 

they  have  gone  around  the  earth.  After  a  wave  goes  through  the 

antipodes  It  returns  toward  the  source.  We  Intend  to  specialise  the 

monochromatic  scalar  wave  equation  In  such  a  way  as  to  eliminate  all 

of  the  waves  which  propagate  from  the  antipodes  back  to  the  source. 

In  other  words  we  will  derive  an  equation  which  Is  first  order  In 

9/36  where  6  Is  the  distance  from  the  source  Instead  of  the 

2  2 

wave  equation  which  contains  9  / 99  .  This  equation  will  be  unable 
to  represent  waves  which  backscatter  from  atmospheric  lnhomgeuelty 
or  waves  which  propagate  more  than  one  half  time  around  the  earth. 
This  equation  does  not  have  periodic  boundary  conditions  but  Is  an 
Initial  value  problem  In  6. 

Before  going  to  spherical  geometry  we  can  more  simply  Illustrate 
the  method  In  cartesian  coordinates.  There  one  wishes  to  restrict 
the  wave  equation  to  waves  with  ray  compnents  In  say  the  plus  z 
direction  rejecting  those  In  the  minus  z  direction.  Let  w  be 
frequency,  c  be  velocity,  and  p  be  pressure.  Then  the  scalar 
wave  equation  In  material  of  slowly  variable  properties  Is  well 
known  to  reduce  to 


A  differential  operator  operating  on  pressure  p  Is  represented  in 
numerical  work  by  a  matrix  operating  on  a  vector  p  where  the  vector 
Is  made  up  of  pressure  sampled  at  various  spatial  positions.  The 
matrix  corresponding  to  a  first  derivative  looks  like,  say 


10 


% 


3  1 

3z  ^  Az 


i 

1  -1 

P(z) 

1  -1  zeros 

p (z+Az) 

1  -1 

p (z+2Az) 

0 

• 

0 

• 

o  ! 

• 

zeros 
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In  mumerical  work  one  o£ten  asks  the  matrix  to  be  as  large  as  possible 

so  that  the  end  effects  may  be  pushed  off  to  infinity.  Then  the 

2  2 

matrix  corresponding  to  3  /3z  is  merely  the  square  of  the  matrix 
corresponding  to  3/3z.  Let  us  rearrange  the  '  tors  in  the  wave 
equation  and  take  a  square  root. 


(2) 


■l 

~  i  2 
Vc 


3x‘ 


17 — 7 

+  l|r=Tr  + 

3x^ 


4 ‘If 


(3) 

(4) 


The  square  root  of  a  differential  operator  may  be  understood  in  terms 

of  the  square  root  of  a  matrix.  A  matrix  A  is  the  square  root  of 

2 

another  matrix  B  if  A  ■  B.  Before  describing  a  systematic  method 
of  finding  the  square  root  operator  we  indicate  why  we  want  it. 
Equation  (4)  is  a  partial  differential  equation  which  is  first  order 
in  z.  Thus  there  is  only  one  independent  solution.  Depending  on 
the  choice  of  sign  of  i  it  is  a  wave  with  ray  componeit  in  the  +z 
direction  or  the  -z  direction.  In  fact  the  meaning  which  we  wish  to 
impart  to  the  positive  square  root  of  the  differential  operator  is 


that  we  will  have  the  same  solutions  as  the  wave  equation  except 
for  those  going  in  a  direction  which  has  a  component  in  the  negative 

z-dlrection.  In  the  customary  approach  to  holography  one  finds  that 

# 

exp(iw(t-r/c))/r  is  the  solution  to  the  scalar  wave  equation.  The 
2  2  1/2 

radius  r  ■  (x  +  z  )  is  defined  by  a  square  root.  The  choice 
of  sign  of  the  root  separates  incoming  from  outgoing  waves.  How 
rather  than  solve  the  wave  equation  and  then  separate  the  solutions 
into  outgoing  and  Incoming  waves  we  are  trying  to  write  an  equation 
which  has  only  the  outgoing  waves  as  a  solution.  If  we  achieve  this 
we  have  converted  a  boundary  value  problem  into  an  initial  value 
problem  and  reduced  the  calculation  accordingly.  To  reassure  the 


reader  that  equation  (4)  is  Indeed  conslstant  with  the  wave  equation 
(3)  one  may  differentiate  (4)  with  respect  to  z 


,  .  |  uf  +  3^  3 

11 1 7  i?* 


(5) 


and  substitute  (4)  into  (5)  getting  the  wave  equation 


The  square  root  of  any  non-defective  matrix  can  always  be  de¬ 
fined  in  an  unambiguous  fashion.  Take  the  matrix  A  to  be  written 
in  diagonal  form  A  -  PAP-1  where  A  is  a  diagonal  matrix  of  eigen' 
values  and  P  and  P  ^  contain  column  and  row  eigenvectors.  Then 
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1/2  1/2  -1 

the  square  root  of  A  is  simply  A  ■  P  A  P  as  may  readily 

1/2 

be  verified  by  substitution.  The  matrix  A  is  created  from  A 

simply  by  taking  the  square  root  of  each  of  the  diagonal  elements 

of  A.  There  Ls  always  a  choice  of  sign  in  a  square  root.  Choosing 

all  plus  signs  amounts  to  choosing  all  waves  in  one  direction. 

Choosing  all  minus  signs  amounts  to  choosing  all  waves  in  the  other 

direction.  Because  of  this  all  one  may  say  that  it  is  conceptually 
22  22  1/2 

clear  how  (to  /c  +  3  /3x  )  has  an  unambiguous  definition  in  terms 
of  matrices  regardless  of  space  variations  in  velocity  c.  This 
definition  is  too  awkward  for  computational  purposes.  A  useful  ap¬ 
proximation  if  the  x-variation  of  velocity  c  has  small  second  de¬ 
rivative  is  given  by  the  binomial  theorem,  the  first  two  terms  of 
which  are  given  by: 

(w2/c2  +  32/3x2)1/2  %  w/c  +  (  c/(  2 to)  )  32/3x2 

Other  more  accurate  approximations  may  be  developed.  A  more  rigor¬ 
ous  derivation  with  computational  details  is  given  by  Claerbout 
(1970). 

Next  we  will  express  these  ideas  on  the  sphere.  The  monochrom¬ 
atic  scalar  wave  equation  is 

(V2  +  to2/c2)  p  ■  0 

On  the  surface  of  a  sphere  the  operators  are 


2  4  2q 
r  sin  0 


sin  6 


30 


H 


sin  0  Is-)  +  — j 

30 '  3<f,2 


co2/c2 


-  0 
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rearranging  and  taking  the  square  root  as  be£ore  we  get 

sin  6  -Tq  "  i((wr  sinG/c)2  +  32/3<f>2)1/2  (6) 

The  square  root  of  the  differential  operator  may  be  defined  as  before. 
One  starts  the  calculation  at  some  0  >  0  by  means  of  an  analytic 
solution  for  the  Lamb  mode.  The  analytic  solution  may  be  determined 
by  the  method  of  Press  and  Harkrider.  Then  by  utilizing  some  approx¬ 
imation  to  (6)  one  could  extend  the  solution  in  the  direction  of  the 
antipodes.  Because  the  velocity  c  may  be  a  realistically  compli¬ 
cated  function  of  space  the  wave  is  expected  to  be  a  mess  as  it  ap¬ 
proaches  the  antipode.  Another  factor  disturbing  perfect  focus  will 
be  the  influence  of  winds.  We  may  avoid  a  long  winded  derivation 
and  get  a  first  order  wind  correction  by  the  following  consideration: 
The  effect  of  a  constant  velocity  wind  on  a  flat  earth  is  exactly 
the  same  as  a  windless  earth  with  a  moving  coordinate  frame.  The 
effect  is  that  an  observer  not  moving  at  the  air  speed  sees  a 
Doppler  shifted  frequency.  Let  w  represent  the  frequency  seen  in 
the  earth  fixed  coordinate  frame  and  ft  represent  the  frequency 
seen  in  the  no-wind  coordinate  frame.  It  is  clear  that  the  exten¬ 
sion  of  (6)  to  include  wind  is 


sin  0  Is-  -  i  ((ftr  sin0/c)2  +  32/3$2)1/2 


(7) 


By  simple  geometrical  considerations  the  two  frequencies  are  related 
by 


where 


ft  ■  w  -  k  •  v 

k  is  the  wave  number  vector  of  a  ray  and  v 


(8) 

is  the  wind 
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velocity  vector.  The  calculation  Indicated  by  (7)  makes  no  provision 
for  definition  of  a  wave  number  vector  1c  however  to  first  order  It 
Is  merely  a  vector  directed  away  from  the  explosion  of  magnitude 
oj/c.  As  the  wind  and  temperature  mess  up  the  expanding  wave  more  and 
more  so  that  the  outgoing  wave  represents  a  wider  and  wider  beam  a  more 
refined  expression  than  (7)  may  be  required.  In  any  event  (7)  is 
inadequate  to  the  task  of  carrying  the  wave  through  the  antipode. 

We  have  not  yet  conceived  of  a  clean  and  neat  procedure  for 
carrying  the  wave  through  the  antipodal  region.  We  have  a  manageable, 
though  messy  solution.  It  depends  upon  splitting  the  collapsing 
wave  Into  several  azimuthal  zones,  transforming  to  cartesian  coord¬ 
inates,  propagating  each  zone  through  the  irregular  focus,  and  then 
transforming  back  to  spherical  coordinates.  We  find  this  unsatis¬ 
fying  in  a  practical  way  but  philosophically  it  encourages  us  to 
beleive  that  there  does  exist  a  good  procedure  of  getting  through 
an  irregular  antipode  and  we  should  be  able  to  work  it  out. 


Figure  4  shows  a  sector  of  a  wave  front  going  through  antipodes  (focus) 
in  the  cartesian  approximation.  The  waves  were  calculated  numerically 
although  an  analytic  solution  was  possible  in  this  case.  Sparsity  of 
data  prevented  a  calculation  based  on  observations. 
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III.  DATA 

A.  Description 

1.  Stanford  data  and  meteorological  influences 

The  records  are  organized  on  a  daily  basis  (from  0:00 
AM  to  12:00  PM).  The  7  records  correspond  I ng  to  the  data  of 
a  week  (In  the  order  of  Sunday  through  Saturday)  are 
collected  in  a  single  sheet.  The  size  reduced  photo  copies 
are  shown  In  Fig. 6.  through  Fig. 13. 

The  records  showed  daily  surface  atmospheric  pressure 
fluctuations  with  relatively  high  frequencies  In  the  after¬ 
noon.  In  the  period  of  6/1  to  10/18/1969,  a  drift  In  the 
mean  pressure  was  observed  dally  from  8:00  AM  to  8:00  PM 
This  drift  was  due  to  afternoon  heating  of  the  Instrument. 
From  November  on.  Dr  Theodore  Madden's  bandpass  filter 
(see  section  111,  B-3)  was  successfully  Incorporated  in  our 
system.  The  high  frequency  fluctuations  and  the  drift  In  the 
mean  pressure  have  been  removed. 

Interesting  events  were: 

1)  It  was  locally  quiet  for  at  least  two  weeks,  but  at  11:30 
PM  6/15/1969,  a  rather  high  amplitude  signal  (about  a  milli¬ 
bar)  arrived.  It  lasted  for  an  hour  and  It  looked  1 1 ke  a 
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dispersed  wave  at  the  beginning,  but  not  so  near  the  end. 

The  largest  period  was  about  10  minutes. 

2)  7/13/1969  2:15-9:00  AM 

A  long  period  signal  came,  its  amplitude  gradually  increased 
and  then  decreased. 

3)  7/141969  3:00  PM 

The  amplitude  suddenly  increased  twice  while  the  frequencies 
were  as  high  as  usual.  It  lasted  for  four  hours.  The  same 
thing  happened  at  3:30  PM  7/24/1969  and  at  4:30  PM  8/3/69. 

The  possible  cause  of  these  events  was  the  'sea  breeze*. 

4)  Starting  at  3:30  AM  9/29/1969,  a  series  of  large  amplitude 
waves  arrived.  The  weather  was  normal  and  quiet  for  about 
three  months.  At  first,  we  misunderstood  i t  as  a  nuclear 
explosion.  However,  a  normal  dispersed  signal  which  lasted 
for  30  minutes  was  found  at  11:05  AM  and  this  was  now  believed 
to  be  a  nuclear  explosion  from  Nop  Nor  (lat:40  degrees  N, 

long:  90  degrees  E)  in  north  west  China.  This  event  was  reported 
in  newspapers.  A  cold  weather  front  with  high  altitude  trough 
existed  some  1500  miles  off  shore. 

5)  A  long  period  (about  8  minutes  in  average)  and  big  amplitude 
(about  a  millibar)  signal  started  at  11:15  PM  10/26/1969  and 
ended  at  5:15  AM  10/27/1969.  We  believe  that  was  the  pass- 
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i  ng  of  the  cold  front. 

6)  At  7:45  12/2/1969,  a  fantastic  large  amplitude  signal 
suddenly  appeared/  It  continued  for  3  hours  and  gradually 
died  out.  There  was  a  cold  front  across  north  western 
America.  The  high  altitude  low  pressure  core  existed  200 
miles  In  the  west.  The  weather  map  and  the  pressure  record 
are  shown  in  Fig. 5. 

7)  1/23/1970  -  2/16/1970 

The  raining  days  have  generally  high  frequencies  and  large 
amplitudes  in  the  pressure  records. 

8)  3/7/1970 

The  solar  eclipse  occurred  today.  The  path  of  the  total  eclipse 
was  from  the  Pacific  to  the  North  Atlantic,  through  Mexico, 
eastern  edges  of  U.S.,  Nova  Scotia,  and  Newfound.  As  a  great 
shadow  swept  the  earth,  the  temperature  dropped  suddenly 
the  effects  on  the  acoustic-gravity  waves  should  be  seen  at 
Stanford  about  4:00  PM  today,  but  from  the  record,  nothing 
special  was  found. 

9)  3/10/1970  7:40  PM  -  3/11/1970  7:30  AM 

There  were  long  period  waves  with  occasional  high  amplitude 
signals.  A  sinusoidal  wave  started  at  6:08  AM  3/11/1970  and 
lasted  for  an  hour. 


a)  The  surface  weather  map  (8:00  PM  12/1/1969). 


b)  The  surface  pressure  record  at  Stanford.  (7:08  PM  12/2/ 
1969  -  1:00  AM  12/3/1969) 


Fig. 5.  The  weather  map  and  the  surface  pressure  record. 

The  possible  cause  of  this  event  Is  by  the  approaching  cold 
front  from  north  west  and  the  high  altitude  Jet  stream  from 
south. 
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III.  A-2  Tel ebarom 


It  was  reported  In  the  newspaper  that  at  8:00  AM  EDT  9/28 
1969,  China  exploded  a  three  megaton  nuclear  device  In  lower 
atmosphere.  The  blast  occurred  In  Lop  Nor  atomic  test  area 
which  we  take  to  be  at  40  deg.  N,  90  deg.  E. 


The  normal  dispersed  signal  arrived  here  at  11:05  AM  9/29 
/ 1 9 6 9  (see  Fig. .14.).  The  longest  period  observed  was  about  8.0 
minutes.  The  event  lasted  for  40  minutes  and  the  periods  were 
gradually  decreasing.  If  the  newspapers  were  right,  the  travel 
time  of  the  acoustic-gravity  waves  due  to  this  explosion  would 
be  6  hours.  The  shortest  path  from  north  west  China  to  Stanford 
is  via  the  Arctic.  In  principle  the  frequency  dispersion  In 
this  record  should  provide  us  with  some  information  of  the 
average  arctic  upper  atmospheric  conditions  in  autumn.  In  fact 
since  there  are  many  unknown  atmospheric  parameters  and  few 
observations  it  is  difficult  to  estimate  an  atmospheric  model 
from  this  data  alone. 


Fig. 14.  The  telebarom  of  a 
1969  Chinese  atmospheric 
nuclear  explosion  (the 
arrow  shows  the  event) 
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IN.  B.  Data  Acqulslons 
1.  The  Playback  Propran  for  LASA  Tapes 

The  LASA  records  of  ml crobaropraph,  weather,  and 
1  onit  period  selsnlc  Information  are  performed  on  an 
incremental  marnetic  tape  recorder.  The  tape  format  Is 
as  follows: 

Bepinninp  of  the  tape  |  TOD  |  Frame  1  I  Frame  2  | . 

I  Frame  120  |  record  pap  |  TOD  |  Frame  1  |  Frame  2.I.. 

I  Frame  120  |  FOF  pap  I  Fnd  of  the  tape 

The  time  of  the  day  (TOD),  conslstinr  of  2  words  and 
end  of  file  (FOF),  are  identical  to  the  LASA  hirh  rate  and 
low  rate  words. 

Fach  frame  contains  data  from  one  sample  time  and  is 
recorded  in  the  follow! np  format: 


Frame  Word 

Site 

Function 

1 

AO 

temp. 

O 

A0 

wind  dir. 

3 

A0 

LTV  baro 

4 

A0 

wind  sp'* 

5 

F3 

LTV  baro 

6 

F  4 

II 

7 

A0 

II 

8 

F4 

II 

Frame  Word  Site  Function 

16  B3  FSSA  bare. 

17  04  " 

18  B4  " 

19  01  " 

20  02  " 

21  B2  " 

22  C3  " 

23  D3  " 
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Frame  Word 

Site 

Function  |  Frame  Word 

Site 

Funct I  on 

10 

FI 

II  j 

25 

D1 

II 

11 

E2 

II  j 

26 

D2 

II 

12 

F2 

II  j 

27 

F3 

LP  Z 

13 

A0 

barometer  j 

28 

F4 

II 

14 

B 1 

ESSA  bar  I 

20 

AO 

II 

15 

A0 

II  | 

30 

E3 

II 

31 

F4 

II  1 

32 

El 

II 

33 

FI 

II  | 

34 

E2 

II 

35 

F  2 

II  | 

Recording  Is  on  7  tracks  at  556  bp! .  Each  record  will 
contain  120  frames  of  data#  each  frame  containing  35  words. 

Thus,  for  3  characters  per  data  word  there  are  12,606  char¬ 
acters/record.  Each  record  Is  then  12,606/556+. 75“23. 45*  Inches 
In  length  and  represents  1  minute  of  real-time  data.  The  amount 
of  data  recorde  I  per  tape  Is  then; 


hours/tape-tape  length  (ft)/117.2  ft/hr. 


recording  for  20  hours  will  require  2344  ft.  of  tape, 
which  Is  about  one  reel. 


The  follow! nr  prorram  ' LASA'  with  JCL  and  comments  Is 
written  for  the  purpose  of  readlnr  the  tapes. 


o  o  o  a  o 
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//LASA  JOB  (J004,315,.9,8,,,,T),'LEE  LU' ,MSGLEVEL»1 

//STEP1  EXEC  FORTHCLG 
//FORT. SYS  I N  DO  * 

THIS  PROGRAM  IS  TO  READ  THE  LASA  TAPES 

REAL*8  DUNCE 
LOG I CAL*1  I NP( 12620 ) 

EQU I  VALENCE ( I  NP(  D/DUNCE) 

I NTEGER*2  X(3, 4202), Y(120, 35), ISIGN, FRAME, INT(8,3,2 ),TOD( 10) 
THE  NUMBER  OF  THE  START  AND  THE  END  OF  THE  RECORDS  TO  BE  READ 
FROM  THE  PROGRAM 

DATA  NRECF/1/, NRECL/2/ 

10  FORMAT( 36Z2 ) 

20  FORMAT( IX/ 1817) 

30  FORMAT(1X/3I1/'j,#2I1/,:,/2I1/,:,/2I1#,.,/I1)  • 

50  FORMAT('l') 

DO  4  KK«NRECF,NRECL 
CALL  NOERRD(  I  NP,  L,  I  ) 

IF(KK.LT.NRECF)  GO  TO  4 
I  F(KK.GT.NRECL)  STOP 
C  3  CHARACTERS  FORM  A  WORD 
C 

DO  5  1-1,4202 
LL-( I -1)*3 

C  IflDK  FOR  THE  TIME  OF  DAY  IN  THE  BEGINNING  OF  THE  TAPF 
C 

DO  5  J»l,3 

5  X( J, I ) -I NP( LL+ J ) 

DO  6  J«l,2 

DO  6  K«l,3 
LL-256 
DO  7  L«l,7 
LL-LL/2 

INT(L,K,J)«X(K,J)/LL 
7  IF(INT(L,K,J).EQ.l)  X(K, J)-X(K. J)-LL 

6  INT(L,K,J)*X(K,J) 

TOD(l)«INT(3,l,2)*2+l  YT(4,1,2) 

TOD(2)-INT(5,l,2)*8+INT(6,l,2)*4+INT(7,l,2)*2+INT(8,l,2) 

TOD( 3)-l NT(3,2,2)*8+INT(4,2,2) *4+  INT(5,2,2)*2+INT(6,2,2) 
TOn(4)«INT(8,2,2)*2+INT<3,3,2) 

TOD(5)»INT(4,3,2)*8+INT(5,3,2)*4+l NT (6,3,2)*2+INT(7,3,2) 

TOD( 6) -I NT(3,1,1)*4+INT(4,1,1)*2+I NT (5,1,1) 


C  LOOK  FOR  THE  DATA  IN  EACH  FRAME 
C 

DO  1  J-1,120 
JBOT ■( J-l)*35+7 
JTOP- JBOT+34 
FRAME-0 

DO  1  K- JBOT,  JTOP 

FRAME-FRAME+1 

ISIGN-1 

I F (X(  1,K)  .GT.3 )  I S I GN  — 1 
I F ( ISIGN.EQ.-l)  GO  TO  2 
Y ( J, FRAME) -X(  1,K) *4096+X( 2 ,K)*64+X( 3 ,K) 

GO  TO  1 

2  Y  (J,  FRAME)  -  (X(  l,K)-60)*4096+X(  2,K)**64+X(  3,K)  -16384 
1  CONTINUE 
WR I TE(6,50) 

WRITE(6,30)(TOD(I),I-1, 10) 

WRI TE(6, 20) (CY(  I, FRAME), FRAME-1, 35, 2),  1  -1,  120) 

4  CONTINUE 
STOP 
END 

/* 

//GO. TWO  DD  DISP-(OLD,KEEP) , VOLUME-SER-LAMA, LABEL-( ,BLP) ,UNIT-TAPE7,  X 
//  DCB-CDEN-1, RECFM-U, BLKS IZE-32000) 

I* 


33 


III.  8-2  The  play  back  procedures  for  Stanford  A/D  tapes 
The  dally  digitized  pressure  fluctuations  are  recorded 
on  the  7  track  tapes  by  our  digital  data  acqulsion  system. 

The  format  of  the  tapes  is: 

Beginning  of  the  tape  I  record  1  |  record  gap  |  record  2  |.... 

I  record  N  |  EOF  gap  |  record  1  |  record  gap  |  record  2| . 

I  EOF  gap  J  end  of  the  tape 

The  data  Is  In  the  form  of  two  bytes  or  12  bits  per  word. 

The  6  bits  In  the  first  byte  together  with  the  first  4  bits  of 
the  second  byte  form  the  numerical  value  of  a  word.  The  last  two 
bits  In  the  second  byte  of  the  word  are  called  the  'u'  bit  and 

the  'm'  bit  respectively.  The  ’u'  bit  has  the  value  1  and  0 

alternately  on  every  other  word,  the  'm'  bit  has  the  value  1  If 

the  data  comes  from  analog  channel  1  and  has  the  value  0 

otherwise.  The  relation  among  the  hexadecimal,  the  decimal  data 
words  and  voltage  of  of  the  output  from  microbarograph  is  shown 
In  the  example  below: 


% 
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Hex 

U_ 

i! 

Decimal 

Vo  1  ts 

0000 

0 

0 

0 

-10.00 

1E2A 

1 

0 

490 

-0.43 

1E2B 

1 

1 

490 

-0.43 

2000 

0 

0 

512 

0.00 

202A 

1 

0 

523 

0.22 

202B 

1 

1 

523 

0.22 

3F3C 

0 

0 

1023 

10.00 

The  end 

of  file 

gap 

(EOF  gap)  can 

be  wr i tten 

the  tape  by  pressing  the  'file  gap'  button  on  our  system#  usual¬ 
ly  we  write  the  end  of  file  gap  once  a  day. 


The  multiplexor  has  8  channels.  So  far  only  the  first  channel 
has  been  used  for  the  output  from  the  microbarograph.  The  number 
of  words  per  record  and  the  speed  of  recording  are  changeable. 
Usually  we  use  256  words  per  record  and  0.05  data  per  second. 


As  an  example  shown  in  Fig,  15.#  good  agreement  has  been 
found  between  the  ml crobarograms  from  the  RUSTRAK  recorder  and 
from  the  plot  of  Stanford  A/D  tape. 

The  computer  program  'tpread*  is  written  for  reading  the 
Stanford  A/D  tapes.  The  second  program  'baro'  is  for  plotting 
the  mi crobarograms  from  the  first  one.  The  computer  programs  with 
proper  JCL  are  shown  in  the  next  few  pages. 
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//TPREAD  JOB  (J004,315,.9,8#/,/T)/ 'LEE  LU' ,MSGLEVEL"1 

//STEP1  EXEC  FORTHCLG 
//FORT. SYS  I N  DO  * 

C  THE  FOLLOWING  PROGRAM  IS  TO  READ  THE  STANFORD  A/D  TAPES 

C  USE  'RUN  HOLD'  TO  SAVE  DATA  ON  DESK 

C 

REAL*8  DUNCE 

LOG  I CAL*1  I NP( 62000 ) 

EQUIVALENCE  ( I  NP(  1),DUNCE) 

I NTEGER*2  INT(512),M(512),U(512),G/H 
C  DEFINE  THE  NUMBER  OF  FILES  ( NF I LES )  TO  BE  READ  IN  THIS 
C  PROGRAM 

NFI LES-5 

WRI TE( 6, 2 )NF I  LES 
2  FORMATC  NFILES-M3) 

IJ  «  1 
NRFC  -  0 
3  CONTINUE 

C  THE  BULIt-IN  SUBROUTINE  NOERRD  IS  USED  HERE  TO  GET  DATA 
C  FROM  EACH  RECORD 

10  CALL  NOERRD  (  I NP,  L,  I  ) 

IF  (  I  .EQ.  1  )  GO  TO  50 
NREC  «  NR EC  ♦  1 
WRITE(6,238) 

238  FORMAT(////) 

I F (  I  . EQ. -1  )  PRINT  9 

9  FORMAT(  '  READING  ERROR  ENCOUNTERED  IN  FOLLOWING  RECORD') 
Ll-L/2 

WRI TE(6,232)NREC, L 

232  FORMATC  RECORD  NO.  ',15,  I  10/  BYTES  LONG') 

C  CHECK  DATA  IN  HEXADECIMAL  FORMAT  BY  USING  THE 
C  FOLLOWING  WRITE  STATEMENT 
C  WRITE(6/22)(INP(K)/K»1,L) 

C  22  FORMAT( IX, 32Z2 ) 

DO  101  K»1,L 
101  INT(K)«INP(K) 

J»1 

K«0 

C  LOOK  FOR  THE  DECIMAL  DATA 
C 

105  DO  102  L«J,256 
K-K+l 
K2«2*K 
H»INT(K2)/4 
G»INT(K2)-H*4 
U(L)-G/2 
M(L)-G-U(L)*2 


n  r>  a  n  o  a  n  r>  n 
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IU-U(L) 

IFU.EQ.  DfiO  TO  111 
I F( ( IU+1 J).NE. 1)  GO  TO  103 
111  JU-U(L) 

102  I  NT <  L ) I NT(K2-1)* 16-512 
GO  TO  106 

103  K-K-l 
DO  104  J-L,256 
K-K+l 

K2AD01-2*K*1 
H-INT(K2ADDl)/4 
G-INT(K2ADD1)-H*4 
U(J)-G/2 
M(J)«G-U(J)*2 
IU-U(J) 

I F ( ( I U+JU) . NE. 1)  GO  TC  105 
JU-U(J) 

104  I  NT ( J  )■!!♦ I NT(K2 ADD 1-1 ) *16-512 
106  CONTINUE 

WRITE(6,235)(INT(K),K-1,256) 

235  FORMAT( IX, 1615) 

check  u  and  m  hits  by  uslnr  the  following  state¬ 
ments 

WR I TE( 6/ 236 ) 

236  FORMAT( /,  '  U:') 

WRITE<6,235>  (U(K),K-1, 256) 

WRITE(6,237) 

237  FORMAT( /, '  Ms*) 

■JRI  TE(6^  235)  <M(K),K-1,256) 
no  TO  3 

50  CONTINUE 
WR I TE(G, 51)  IJ 
NR EC  ■  0 

51  FORMATCOEND  OF  FILE  NO.*, 14,/,  '1',//) 

I F (  IJ  .GF.  NFILES  )  STOP 
IJ  -  IJ  ♦  1 

C  PROGRAM  AROUND  A  BUG  IN  NOERRD  BY  MAKING  AN  EXTRA  CALL  MOW. 

CALL  NOERRD(  I  NP,  L,  I  ) 

GO  TO  10 
END 

/* 

//GO. TWO  DD  DISP=(OLD,KEFP),VOLUME-SER»GOOO,LABFL=(,BLP),UNIT»TAPE7,  X 
//  DCB=(RFCFM=U,BLKS1ZE=32000) 

/* 


oooor>oo  oo  oo  ooo 


V 
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//BARO  JOB  (J004,315,  .9,2,,,,T),  ’LEE  LU * ,MSGLEVFL«1 
//  EXPO  PLOTCLG 
//FORT.SYSI N  DD  * 

THE  FOLLOWING  PROGRAM  IS  TO  PLOT  THE  OUTPUT  FROM  LASA 
OR  STANFORD  A/D  TAPES 

I NTEGER*2  I NT( 1536 ) 

REAL*4  R I  NT C 1 5 36 ) 

READ  DATA  FROM  OUTPUT  OF  THE  TAPE  AND  CONVERT  TO  REAL  NUMBERS 

READ( 5/1)  (INT(I),  1=1,  1536) 

WR I T  E  ( 6  ,  1 )  (INT(I),  1=1,  1536) 

1  FORMAT( 1615) 

DO  2  1=1,1536 
J  =  l  NT(  I  ) 

2  RINTCI )«FL0AT(J) 

WRITE  THE  TITLE  AND  PLOT  THE  TIME  SERIES 
CALL  STRTPl(lO) 

CALL  SYMBL1C3. ,7., .3, 'MICRO BARO.  6  RECORDS  START  AT  10:50  PM  5/12 
* / 6  9 ' , 46 ) 

CALL  PLOTl( 0. ,2.5,23) 

CALL  TPLOTC  1,  1536, Rl NT,  1. , 8 . 5, 1. , 1, 1. ) 

CALL  ENDP1 
STOP 
END 

SUBROUTI NE  TPLOT( LI MES, L, X, F I LL, XI NCH, Yl NCH, MEANS, PTDELT) 
FILL-OVERLAP  FACTOR:  1  IS  NO  OVERLAP, 10**6  FOR  COMMON  AXIS 
XltyCH-LENGTH  IN  INCHES 
YINCH-WIDTH  IN  INCHES 

IMEANS | =1  FOR  MEAN  LINES  AND  MEANS. LT.O  FOR  VERTICAL  AXIS 
PTDELT-SEPARATION  FROM  PRECEEDING  PLOT 
IF  FILL  IS  NEGATIVE  THFN  PIPS 

MUST  PLACE  CALL  STRTPl(lO)  AND  END PI  BEFORE  AND  AFTER  CALLS  RESPEC 
DIMENSION  X  C  L ) 

OPT-S I GN( 1. , F I LL) 

M»  I S  I  GN(  1,MFANS ) 

MEANS-1 ABS (MEANS) 

FILL-ABS(FILL) 

DTI  N-AMAXK  .01,  (X I  NCH*LI  NES ) / L) 

SIDES-1. 

B-0. 

DO  5  I =1, L 

IF(X(I ).LT.0)SIDES=2. 

5  I F  C  ABS ( X ( I ) ) .GT. B)B*ABS( X( I ) ) 

DY=F  I  LL*Y  I  NCH/  (F I  LL+LI  NES -D/SIDES 

DM=S I DES*DY/F  I  LL 

CALL  PLOT 1( PTDELT,  5. , 23) 

K-L/LI NES 


SP-YINCH/2.-DY 
DO  3  I -1,  LINES 
CALL  PLOTHO. ,SP,23) 

IF(M.NE.-l)GO  TO  2 
CALL  PLOTHO.  ,DY,  2 ) 

X1*0. 

IF(SIDES.EQ.2)X1— DY 
CALL  PLOTHO.,  XI,  2) 

CALL  PLOTHO. ,0. ,3) 

2  CONTINUE 
SP— DM 

YaX( 1+K*( I -1) )*DY/B 
Z»0. 

CALL  PLOTHZ,  Y,  3) 
if(opt)8,8,9 

8  CALL  PLOTHZ/  0.  /  2 ) 

9  CONTINUE 
DO  4  J«2,K 
Y«X(J+K*(I-1))*DY/B 
Z-DTIN+Z 

I F (OPT) 7/  7,6 
7  CALL  PLOTHZ, 0., 2) 

CALL  PLOTHZ,  Y,  2 ) 

CALL  PLOTHZ, 0., 2) 

GO  TO  4 

6  CALL  PLOTHZ, Y, 2) 

4  CONTINUE 

Z-Z+DTIN 

IF (MEANS. NE. l)GO  TO  3 
CALL  PLOTHZ, 0., 3) 

CALL  PLOTHO. ,0., 2) 

3  CALL  PLOTHO. ,0., 3) 

IFCSIDES. EQ.l. )DY-0. 

DY  —  l.*(DY+5.-YINCH/2.) 

XI MCM “X I NCH-DTI N 
CALL  PLOTHX I  NCH,  DY,  23) 

RFTIJRN 

END 

/* 

//GO. SYS  IN  DD  * 

(THE  DATA  FROM  TAPES  IN  THE  FORMAT  OF  1615) 
/* 


Fir.  15.  :  M i crobarorram  (start  from  10:4b  PM  12/4  to  1:35  AM  12/5/59) 

a)  Tho  record  from  RU5TRAK  recorder,  (1  inch/hour). 

b)  Tho  plot  from  Stanford  A/D  tape.  (2  records  or  512  words 


of  data  in  tho  lon.cth  of  G  inches,  the  speed  of  the 
record! nr  is  20  seconds/data  word). 


V 


* 
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III.  B-3  The  Ml crobaropraph  Filter 

The  Information  of  the  band-pass  filter  was  provided  by 
Dr.  Theodore  Madden.  The  hlph-pass  time  constant  Is  500  sec 
equivalent  to  45  min.  cut  off  which  Is  somewhat  longer  than 
the  leak  cutoff.  The  low  pass  Is  not  critical  due  to  the 
natural  fall  off  of  the  pressure  fluctuations  with  frequency 
for  periods  less  than  the  Brunt  Valsala  period. 

The  circuit  diagram  of  the  filter  Is  following: 


50  M 


Fig. 16.  The  circuit  diagram  of  the  bandpass  filter 
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The  'A'  In  Fig. 16.  Is  a  low  current  offset  and  high  Impedance 
FET  or  parametric  amplifier.  We  used  the  FET  operational 
amplifier  Model  141A  by  ANALOG  DEVICES  INC.,  because  of  the  low 
cost  and  the  low  Input  bias  current. 

The  gain  and  the  pass  band  are  determined  by  the 
character  I  sties  of  the  equation: 

V(out)  —  ((R+Rg)/R)*(Z2/Zl)*V(ln) 

Where,  the  gal n«(R+Rg)/R,  and  In  our  case 

Z2/Zl«500.*J*w/((l.+5.*J*w)*(l.+500.*J*w))  • 

The  pass  band  Is  calculated  from  two  poles  In  Z2/Z1,  le.  the 
angular  frequency  Is  In  between  0.2  and  0.002  cps  or  the  period 
Is  within  the  passband  of  0.5  to  52.3  minutes.  The  characterlc 
curve  of  I Z2/Z 1 |  Is  shown  In  Fig. 17.  An  example  of  the  filtered 
and  the  unflltered  pressure  records  are  compared  In  Fig. 18.  The 
high  frequency  fluctuations  00.2  cps)  have  been  successfully 
removed. 


42 


1  .  'j  +-  --  - 


1.0  + 


|Z2/Z1| 


0. 


r, 


****** 


*l< - - --->1 

tho  pnssband 

*  * 


o.n  +-• 


io"4  io~3  io~2  io”1  io( 

the  angular  frequency  (rps) 


10 


Fir.  17.  Thr*  character  i  s  t  i  c  curve  of  |Z2/Z1| 
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Fig. 18.  a)  The  unfit tered  microbarogram. 
b)  The  filtered  m 1 crobarogram. 
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